options(scipen=999)

rm(list = ls())
library(data.table)
library(rdrobust)
library(lmtest)
library(multiwayvcov)
library(stargazer)

setwd(dir = '.') # Please insert folder where files are stored
load('ver.RData')

rural.data <- ver[ruralPop > median(ver$ruralPop,na.rm=TRUE)]
  bw.r <- rdbwselect(rural.data$f_elected, rural.data$marginal, bwselect  = 'msetwo', masspoints = 'off')
rural.data.a <-rural.data[abs(marginal) < bw.r$bws[[1]] & elected == 0] # 0.004618702
rural.data.b <-rural.data[abs(marginal) < bw.r$bws[[2]] & elected == 1] # 0.01130364
rural.data <- rbind(rural.data.a, rural.data.b)


weights <- kernelwts(rural.data$marginal,center = 0, bw = max(bw.r$bws[[2]],bw.r$bws[[1]]), kernel="triangular")
rural <- (lm(f_elected ~ I(inflation*elected) + inflation + elected*marginal + I(inflation*marginal)+as.factor(year),data=rural.data))
  vcov_both <- cluster.vcov(rural, cbind(rural.data$muni_code))
  rural.coeftest <- coeftest(rural, vcov_both)

urban.data <- ver[ruralPop < median(ver$ruralPop,na.rm=TRUE)]
bw.u <- rdbwselect(urban.data$f_elected, urban.data$marginal, bwselect  = 'msetwo', masspoints = 'off')
urban.data.a <-urban.data[abs(marginal) < bw.u$bws[[1]] & elected == 0] # 0.002394197
urban.data.b <-urban.data[abs(marginal) < bw.u$bws[[2]] & elected == 1] # 0.01017874
urban.data <- rbind(urban.data.a, urban.data.b)

weights <- kernelwts(urban.data$marginal,center = 0, bw = max(bw.u$bws[[2]],bw.u$bws[[1]]), kernel="triangular")  
urban <- (lm(f_elected ~ I(inflation*elected)+inflation+elected*marginal+I(inflation*marginal)+as.factor(year),data=urban.data, weights = weights))
  vcov_both <- cluster.vcov(urban, cbind(urban.data$muni_code,urban.data$year))
  urban.coeftest <- coeftest(urban, vcov_both)

stargazer(rural.coeftest,urban.coeftest,stars='default',
          digits=3,
          type='text',
          keep.stat = 'n',
          add.lines = list(c("n",rural$df.residual, urban$df.residual)),
          column.labels   = c("Rural","Urban"),
          omit = c('year','marginal'),
          dep.var.caption = 'Probability Win Next Election', dep.var.labels = '',
          label = 'vereadores',
          covariate.labels = c('$\\Pi_{ki_{t+1}} \\times Incumbency_{ki_{t}}$','$\\Pi_{ki_{t+1}}$','$Incumbency_{ki_{t}}$','$Margin$','$\\Pi_{ki_{t+1}} \\times Margin$','$Incumbency_{ki_{t}} \\times Margin$','Constant')
          )

